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Abstract 

The triiodide ion is an example of a system where symmetry-breaking 
may be induced by a solvent. The Landau free energy is expected to have 
a similar form to that for the mean field Ising model, but with solvent 
strength rather than temperature as the control parameter determining 
whether there is symmetry breaking. In order to examine the possibility 
of critical phenomena we have studied the properties of the ion in solvents 
based on a model for water with charges scaled by a factor A. As A is 
increased the system changes from one with no symmetry breaking to 
one with strong symmetry breaking. The behavior of various quantities, 
including the Shannon entropy, as a function of A show only weak maxima 
near the critical value of A = Ac. We examine the behavior of a simple 
model and show that divergences would only be expected in the limit of 
low temperatures, and the essential difference between the solvent-induced 
symmetry breaking and the mean field Ising model is that in the latter 
the observed quantity is an average over many spins, while in the former 
observations are made on individual molecules. 



1 Introduction 

In an earlier paper we showed that symmetry breaking could be induced in 
the triiodide ion by varying the solvent 1 . Experiments and simulations 
m El El C] suggest that protic solvents which can form hydrogen bonds with a 
negative ion cause symmetry breaking of the ion, so that the charge becomes 
concentrated at one end of the ion and the corresponding bond elongates. We 
suggested that one could draw an analogy between the mean field Ising model 
with free energy 

F = a{T - To)7f -f Cry4 (1) 
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and solvent induced symmetry breaking with 

F^a{Ro-RW + Crj^. (2) 

In these expressions F is the Landau free energy per molecule (or per spin), rj 
is an order parameter and a and C are coefficients The Ising model has 
a critical point when the temperature T reaches Tq; below this temperature 
the symmetry is broken to give domains with non-zero order parameters, while 
above this temperature the order parameter is equal to zero. Various critical 
phenomena, such as vanishing of the susceptibility and diverging fluctuations in 
the order parameter, are observed as the critical point is approached. As these 
mean-field critical phenomena result from the form of the equation for the free 
energy, it is interesting to see whether there are corresponding phenomena in the 
solvent-induced symmetry breaking. In this paper we describe an investigation 
of this point. In order to vary the solvent strength in a systematic way we used 
a series of modified waters as solvents. The models are based on the standard 
spc/e model with a Lennard- Jones center on the oxygen atom and charges on the 
atomic sites. The charges are scaled by a factor A which varies from unity, giving 
the standard spc/e model which we know causes strong symmetry-breaking, to 
zero, giving a pure Lennard- Jones solvent which our earlier work shows does 
not cause symmetry breaking. Thus varying A provides a method of tuning the 
solvent strength R through the critical value Rq. 

2 Theory and computational details 
2.1 Valence bond model for I3 

The model for triiodide ion is the same as used in our previous work mSlEIIHlini- 
The electronic structure of the ion is described using a semi-empirical valence 
bond model based on diatomics in molecules method {.iPj with additional terms 
due to the fact that the species is charged . Provided the ion is constrained 
to be linear, its ground state is described by a 3 x 3 Hamiltonian matrix, whose 
matrix elements depend on the instantaneous values of the bond lengths and the 
instantaneous external electrostatic potential due to the solvent. Full details are 
given in the appendix to reference [H] . The solvent molecules are rigid three-site 
models based on spc/e water J2], with partial charges on all three atomic sites 
and Lennard- Jones interactions on the oxygen site. The total energy of the 
system can be written as 

E{{r^}) = E '^o.cof.H^p + J2 -^^^ + E + E (3) 

a/3 ij ^ ^ ^ jk jm 

where the first term is the quantum mechanical energy of the ground state, 
with Haf} being the Hamiltonian matrix element between basis states a and /3 
and cqq the coefficient of basis state a in the ground state. The second and 
third terms are sums of the electrostatic and Lennard- Jones interactions over 
all pairs of solvent sites j and k, and the last term is the sum of Lennard- Jones 
interactions between solvent sites and iodine atomic sites m. 

Quantum mechanical forces on both iodine and solvent sites were calculated 
using the Hellman-Feynman theorem, while the forces arising from the three 
classical terms in equation were calculated in the usual way within the 
molecular dynamics program. The ion was constrained to be linear throughout. 
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2.2 Modified water models 

The solvent molecules were based on the standard spc/e model for a water 
molecule jl2| with a Lennard- Jones center on the oxygen and charges on the 
atomic sites. In spc/e water hydrogen bonds and all other orientational corre- 
lations are due solely to electrostatic interactions between molecules. Thus the 
hydrogen-bond strength can be controlled by scaling the charges. In this study, 
nine solvent models were used with the charges were scaled by a factor A which 
varied from 1.125, giving a super strong water, through 1, giving the standard 
spc/e model, to zero, giving a pure Lennard- Jones solvent. 

The Lennard-Jones potential between the iodine and oxygen sites was the 
same as used in our previous work with water 01 1^ . The Lennard-Jones 
parameters and the partial charges for sites are given in Tabled 

2.3 Simulation details 

Molecular dynamics simulations were carried out using a version of the dlpoly 
program which was modified to include the construction and diagonalisa- 
tion of the Hamiltonian matrix and the calculation of the Hellman-Feynman 
forces. The simulation cell contained one triiodide ion and 509 modified water 
molecules in a cubic box with an edge-length of 24. SA. The density was the 
same for all simulations. Simulations were carried out at 300K with a time 
step of 1 fs. For each value of the charge scaling constant. A, the system was 
thoroughly equilibrated before collecting data for 1 ns. Our earlier work showed 
that two order parameters were necessary to describe the symmetry breaking. 
The order parameters chosen are the molecular dipole moment fi (relative to 
the center of mass) and the normal coordinate for the antisymmetric stretch 
C = (^12 — b23)/VQ- The former measures the extent of electronic distortion, 
while the latter measures the geometrical distortion. The instantaneous values 
of these parameters were determined at each time step, and averages, mean 
square fluctuations and probability histograms were constructed. 

2.4 Shannon's information entropy 

From the two-dimensional histograms for the probability distributions of the 
two order parameters, the Shannon information entropy function -ff (C, ^J■) was 
constructed for each bin 

i7(C,M)--p(C,Ai)lnp(C,M), (4) 

where /x) is the probability of being in that bin, with ^ X]p(Ci — 1- The 
Shannon entropy of the system ^3] is then 

We show in the appendix that the limit of this sum as the bin sizes tend to 
zero is a property of the system, rather than of our information about it, and 
is the entropy associated with the spread of order parameter values. In our 
calculations the bin sizes used were (5/i — 3D and 5C — O.OAA. The values of the 
bin sizes affect the zero of entropy, but, as they were kept constant for all the 
simulations, entropy differences between runs are real although absolute values 
are arbitrary. 
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3 Results 



Figure n shows the distribution of the Shannon entropy function H{(,fi) for 
selected values of A. When A is small there is a single maximum while when 
lambda is large enough two maxima are seen. Figure [21 shows cross sections 
through this H surface for all the runs. The direction of cross section is dif- 
ferent for each A and is either chosen to go through the two maxima, or, when 
there is only a single maximum, it is chosen to go through the direction of 
minimum curvature. These figures show that there is symmetry breaking in 
the curve marked 4 (A = 0.5) and there is no symmetry breaking in the curve 
marked 3 (A = 0.375) Thus the critical value of A, Ac, lies between these values, 
that is somewhat below A = 0.5. These results confirm that varying the scaling 
parameter A induces symmetry-breaking. 

There is no a-priori reason to identify A as opposed to some function of A 
with the solvent strength in the free energy equation (0) . However figure 13 
shows the solvent-triiodide interaction energy (lowest curve). The interaction 
energy decreases smoothly with the scaling factor A and is approximately linear 
in the critical region. This figure also shows (upper curve) the energetic cost of 
polarising the ion. 

The inverse susceptibilities for the response to an external electric field 

= kT/ <^i^> (6) 

is plotted in figure 01 It can be seen that there is a monotonic decrease in the 
inverse susceptibility as a function of the scaling factor A, and that does not 
tend to zero near the critical value A = Ac. In a second-order phase transition, 
however, inverse susceptibilities do tend to zero at the critical point. 

The Shannon entropy, hs, would also be expected to show critical behavior 
at the critical point. Figure shows that in our system this quantity has a 
weak maximum at a value of A which is slightly greater that the critical value 
Ac, rather than diverging at Ac. 

4 Discussion 

As the two figures O and 01 do not show the expected critical behavior near 
the critical value of the scaling constant, A = Ac, we conclude that there is 
a difference between our system and the mean field Ising model. In order to 
elucidate this difference we study the properties of a simple model. 

4.1 A simple model 

Let us consider a simple model for symmetry breaking with a single order pa- 
rameter r] and a Landau free energy per molecule given by 

F = a{Ro - R)f + Cfj\ (7) 

This can be rewritten (by rescaling F = CF/{aRo)^ and 77^ = Cf]'^ / {2aRo)) as 

F=(l-p)r;V2 + W4, (8) 

where p = R/Rq is a measure of the relative strength of the solvent interaction. 
The critical point where symmetry breaking occurs is p = 1. In order to answer 
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the question as to when one should observe critical phenomena such as diverging 
fluctuations, we examine the properties of the probability distribution of the 
order parameter 

p(r,) = exp[-/3F]/Z - exp[-/3((l - p)ff /2 + ii^/A)]/Z (9) 
as a function of the parameter /3. In this expression 

[ eM-m~ pW/'i+vVmv- (10) 



Figure El shows the values of the entropy h = — J plnpdr] as a function of 
p for a number of values of /3. It can be seen that the Shannon entropy has a 
maximum as a function of p. As /3 increases, this maximum gets sharper and 
shifts downwards towards p = 1- 

For a second order phase transition, the inverse of the susceptibility x, which 
is given by = ^/Piv'^), tends to zero at the critical point where the phase 
transition occurs. The upper part of figure [7| shows the values of the inverse of 
the susceptibility of the model system as a function of p for a range of values 
of the parameter /3. Again we observe that if f3 is large enough the inverse 
susceptibility decreases linearly towards a value of zero at p < 1. A simple 
calculation of the susceptibility above the critical value of p is misleading as the 
probabilities of positive and negative values of the order parameter are equal 
and the mean value is always 0. However if a small biassing field b is added 
giving the free energy 

= -6?7 + (1 - p)?7V2 + ??V4, (11) 

then the lower part of figure shows that the inverse susceptibility increases 
with p above the critical value p = 1. 

The results at large (3 are very similar to those obtained from a mean field 
model of a second order phase transition with Landau free energy 

F ^ -br] + a{To-T)'n^/2 + Cjf/A. (12) 

In the standard treatment of this model j2], one is only concerned with the 
minima of the function. One finds that above the critical temperature Tq the 
minimum of F is at rjrnin = in the absence of an external field, while below Tq 
the minima are at rjmin = ±(a(Io — T)/CY^^. The susceptibility is calculated 
from the change of the position of rjmin with the strength of the external field, 
X = drjrain/db. This gives 

X-^=a{T-T^) T>To 

= 2a{To -T) T< Tq. (13) 

It will be noted that the standard treatment does not include fiuctuations in 
the order parameter, while the expression that we have used (equation inj gives 
the full range of possible values. The reason that one can ignore fluctuations 
in the treatment of phase transitions is that the expression for the free energy 
given in equation (|12|l is the free energy per unit cell or per spin. The observed 
f] is the average over all the unit cells or all the spins and, as one observes the 
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average over a large number N of unit cells or over N spins, the probability of 
observing a given value of 77 is given by 



In the limit A'^ — )■ 00 it is indeed only the minima that are observed. In the 
solvation-induced symmetry breaking situation we observe individual molecules 
and the relevant probability is given by the similar equation IHl but with the 
difference that /3 = {AaRoY /{CkT) rather than N/kT. The fact that we only 
see significant maxima in the Shannon entropy near the critical value of p when 
(3 is large is consistent with the fact that true critical phenomena only occur in 
the limit of — > 00 and depend on the observation of an order parameter which 
is an average over many replicated systems. While solvent-induced symmetry 
breaking will never show true divergences, there will be maxima in the Shannon 
entropy if the value of the parameter {aRo)^/C is large enough compared to 
kT. Comparing the model results with those observed for the triiodide ion, 
we estimate that the order of magnitude of {aRo)"^ /CkT is about 100. At this 
value the is a weak maximum in the entropy at a higher solvent strength than 
the critical one and the inverse susceptibility decreases smoothly through the 
critical point. 

5 Conclusion 

In this paper we have examined the solvent-induced symmetry breaking induced 
by water and modified water. The use of scaled charges in the solvent models 
allows us to vary the solvent strength which induces the symmetry breaking 
continuously and to determine whether there are any phenomena analogous to 
critical phenomena. In this particular system the transition from no symmetry 
breaking to symmetry breaking is weak and there is only a small maximum in 
the Shannon entropy. By comparing the results to a simple model we see that if 
(3 — {aRo)^ /CkT « 100 the model system shows a rather similar behavior. In 
order to see significant critical behavior this parameter would need to be larger 
by a factor of 10 or more. 

Although this example of solvent-induced symmetry breaking does not show 
critical effects, there may be other situations which do. A system which was 
less polarisable, would have a larger value of ai?o, but would need a strong 
interaction with the solvent for symmetry-breaking to occur at all. 
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Appendix 

Standard methods of statistical mechanics shows that in contact with a heat 
bath the probability p(r)dr of observing the 3N system coordinates r = {r^} 
with values between r and r -I- dr is given by 



p{v) = exph(A^/fcr)(a(To - T)7jy2 + Crjy4)]/Z. 



(14) 



p(r) = exp[-(3H{r)]/Z, 



(15) 
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where H is the potential energy. Z is the configurational integral for the com- 
plete system 

Z = j exph/3i7(r)] dr = exp[-/3Feonf], (16) 

which defines -Fconf , the classical configurational part of the Helmholtz free en- 
ergy of the system. 

The total Helmholtz free energy can be written 

F = -kT\nZ-kT In A = F.^ni + F^om (17) 

where Fmom = —kTlnA is an additional ideal term due to the momentum and 
indistinguishability, 

A = J](27rm,fcr//i2)3/2/7vj (18) 

i 

where the product is over different types of atom, i. 

Let us define the Landau free energy F^ (jj) as the free energy of the system 
when an order parameters 77 (r) is constrained to a fixed value rjQ. Thus 



(19) 



From equation (|15(l we see that the probability p{r])d'r] of observing the value 
of the order parameter between 77 and rj + dt] in the unconstrained system is 

p{f]) dr] = exp[-/3Fi (77)] df]/Z (20) 

with 

Z = J eM-f3FL{v)]dv. (21) 
Now let us consider the quantity 

p{ri)liip(r]) drj. (22) 

h is obviously a property of the system (molecule plus bath) rather a property 
of our knowledge about it. In the words of Denbigh and DenbighjTl] it is an 
objective property rather than a subjective property. As p(r]) is a density rather 
than a dimensionless property, Z has the dimensions of the order parameter and 
so the zero of h depends on the units of rj. 
Using the expressions above 



or 



p{r,)FL{ri)dr] - F, 



conf 



(23) 



/l = /3[(Fi(77)) -feonf], (24) 

the difference of the average Landau free energy and the configurational free 
energy of the unconstrained system. However as, F = U — TS and the average 
energy is the same whether calculated by averaging over the constrained sys- 
tems or over the unconstrained system, h is the difference of the entropy of a 
constrained system averaged over all constrained systems and the total entropy, 
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h=[{SLiv))-Scon{]/kB. (25) 

Here 5*^ is the entropy of a constrained system with fixed 77, and the angular 
brackets denote an average over all constrained systems weighted with their 
probabilities in the unconstrained system. Thus h can be interpreted as the 
entropy associated with the order parameter distribution. 

The Shannon information entropy is defined in terms of the information that 
one has about the system. For example if a histogram of probabilities of r] is 
constructed with bin widths S and if one calculates the Shannon entropy using 
natural logarithms rather than logarithms to base 2, 

hs = -^Pilnpi, (26) 

then 

hs ~ J plnpdr] — ln{S) 

= h-\n(S). (27) 

Thus the Shannon entropy, hs, which is a property of the information col- 
lected about the system, provides an approximation to the system property h. 
The latter is a measure of the spread of the order parameter. This treatment 
may readily be extended to more than one order parameter. 
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Table 1: Intermolecular Lennard- Jones site-site parameters for water. Cross 
terms in the Lennard-Jones interactions were calculated using the Lorentz- 
Bertholet rules. 



atoms i 


Cii /kJ mol ^ 




q(i)/e 





0.6502 


3.169 


-0.8476A 


H 








0.4238A 


I 


0.4184 


5.167 


varies 
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Figure 1: Three-dimensional plots of Shannon's information entropy function 
H{(^, ji) of at 300 K as a function of the antisymmetric vibrational normal 
coordinate and the dipolc moment /i for different modified water solvents. 
The distributions correspond to charge scaling factors A = 0,0.375,0.75, 1.125 
(from top to bottom, and left to right). Note the gradual changes from a single 
peak (symmetry preserving) to double peaks (symmetry breaking) and also the 
changes of distortions around each peak. 
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Figure 2: Cross sections through the Shannon entropy function surfaces for 
different A values. The curves are labelled with t = 8A. Note that symmetry 
breaking first occurs when £ = 4 (A = 0.5). 




_io r 1 
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Figure 3: Solvent-solute energetics. The upper curve (triangles) shows the 
change in the average internal energy of the molecule relative to the gas phase 
which is a measure of the cost of polarising the molecule. The lowest curve 
(squares) shows the average solute-solvent interaction energy and the middle 
curve (circles) shows the sum of these two energies, as a function of A. 
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Figure 4: Inverse susceptibility of the triiodide ion in solution as a function of 
A. 
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Figure 5: Total Shannon entropy /i as a function of A for a triiodide ion in 
solution. Note that there is a weak maximum at about A = 0.6 although the 
critical value of A for symmetry breaking lies between A = 0.375 and 0.5. 
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Figure 6: Variationof the entropy associated with the order parameter, h, with 

relative solvent strength p for different vlaues of the parameter /?. Note that 
as P increases the maximum gets sharper and moves closer to the critical value 
p=l. 
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Figure 7: Inverse susceptibility of the model as a function of the solvent strength. 
Above: rcults for various vahics of the papcramctcr (3. Below: Rcsiilts for a 
high value of /3 with and without a small biassing field. Note that with the 
biassing field the model has a minimum inverse susceptibilty near the critical 
value p = 1 . In the absence of the biassing field there is a linear decrease in the 
inverse susceptibilty between p = and p = I. 
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